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Evolution of bimodal accretion flows 
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ABSTRACT 

We consider time-dependent models for rotating accretion flows onto black holes, 
where a transition takes place from an outer cooling-dominated disc to a radiatively 
inefficient flow in the inner region. In order to allow for a transition of this type we solve 
the energy equation, both, for the gas and for the radiation field, including a radiative 
cooling flux and a turbulent convective heat flux directed along the negative entropy 
gradient. The transition region is found to be highly variable, and a corresponding 
variation is expected for the associated total luminosity. In particular, rapid oscillations 
of the transition radius are present on a timescale comparable with the local Keplerian 
rotation time. These oscillations are accompanied by a quasi-periodic modulation of 
the mass flux at the outer edge of the advection-dominated accretion flow (ADAF). 
We speculate about the relevance for the high-frequency QPO phenomenon. 
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1 INTRODUCTION 

The stan dard model for an optically thick, Keplerian accre- 
tion disc iShakura fc Sunvaevll97a) explains the continuous 
emission observed from a variety of accretion sources both 
on the galactic and extragalactic scale. This so-called SSD 
model cannot explain, however, the hard spectral component 
observed in most compact accreting systems, e.g., in X-ray 
binaries or in AGNs. This component is generally fitted by 
a power-law and ascribed to a thin, radiatively inefficient 
plasma. Such conditions are met in an advection dominated 
accretion flow (ADAF, for a review see lNaravan et aT1ll99ct) 
where the heat, that is liberated locally by a disc annulus 
due to viscous dissipation, is advected with the radial bulk 
flow and eventually disappears into the central black hole, 
without being emitted towards the observer. 

The presence of reflection components and the Fe Kg 
line o bserved in X-ray spectra of BHXB's febisawa et alJ 
f 996) strongly suggest that cold gas from a disc and hot gas 
must coexist in a small volume. A possible configuration is 
the radial transition from an outer SSD to an inner ADAF. 
A bimodal disc model along this line, with the transition ra- 
dius as a free parameter and with the ADAF and SSD parts 
treated as separate en tities, allows to fit the spectra of per- 
sistent X-ray sources jNaravan et al.1119961 iManmoto et alJ 
Il997l: lEsin et al]l200ll) . Furthermore, this model does also 
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explain the different spectral states observed in some soft 
X-ray t ransients by vary ing the location of the transition 
radius fesin et alJ If 998). A number of these sources do 
also show high-frequency quas i-periodic oscillations (QPO) 
in their fourier-power spectra jNowak fc Lehrllf 9981) . These 
examples demonstrate the necessity for a time-dependent 
treatment. 



In a consistent theory, the disc transition radius is 
obtained naturally as a function of flow parameters and 
boundary conditi ons. For the stationary l imit, this issue 
was addressed by IManmoto fc Katd (|2000) . If the transi- 
tion region between the inner ADAF and the outer SSD 
is sufficiently narrow (i.e., SR/R <C 1), specific connection 
conditions can be derived which in turn significantly con- 
strain the parameter space available for global solutions. 
The most restrictive condition is that the ADAF must sup- 
port a minimum outward energy flux into the SSD, which 
is then liberated as radiation. The nature of this energy 
flux is still a matter of debate, in particular within th e 
one- dim e nsion al app roximation ll Abramowicz et all l2000f) . 
iHonmal £996) and IManmoto fc Katd i2000t) specifically 
considered energy transport by turbulent convection, and 
conclude that the transition radius for stationary flows is 
then constrained to a radial domain with boundaries de- 
pending on the model parameters. 
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2 THE DISC MODEL 

We investigate time-dependent transition models by solving 
the full Navier-Stokes equations for a vertically integrated 
axisymmetric disc. In the energy equation for the gas we in- 
clude dissipation of a turbulent convective energy flux, and 
in addition we solve a radiative transfer equation in the dif- 
fusion approximation. Our dis c model thus generalises the 
stationary transition model of iManmoto fc Katol i200(tl to 
the non- stationary case. 

Assuming vertical hydrostatic equilibrium, the vertical 
pressure scale height is H — Cg/fiic, with c s the isothermal 
sound speed and Qk the Keplerian angular velocity associ- 
ated with the pseudo-Newtonian gravitational potential for 
the space-time of a non-r otating black hole with mass M 
jPaczvnskv fc Wiital fl980h . As a result of the vertical av- 
eraging procedure fe.g.. iNaravan fc Yilll995al) . the thermo- 
dynamic variables which appear in the model equations are 
all vertically integrated quantities, such as surface density 
S = 2Hp, total pressure P, internal energy density E and 
radiation energy density U . 

The continuity equation and the momentum equations 
in radial and azimuthal direction are 

<9 t E + V-£w = (1) 

dtXvR + V-XvRV = -<9i?P + E(fi 2 -C&)R+(V-t)r (2) 

d t U + V-E/w = (V-r)* R, (3) 

where R is the cylindrical radius and I is the specific angu- 
lar momentum. Viscous transport is modeled as a Navier- 
Stokes stress, i.e., the vertically integrated stress tensor 
t — 2pa, with a the shear tensor and p the vertically inte- 
grated dynamical viscosity . We adopt the a- viscosity law of 
IShakura fc Sunvaevl il973f) as a model for anomalous turbu- 
lent viscosity, v = ac s H ', which results in p = T,v = T,ac B H. 

The thermal disc structure is governed by the magni- 
tude of the various heating and cooling terms which con- 
tribute to the energy equation for both the gas and radiation 
field. The gas energy equation is written as 

d t E + V-EV + P g V-v = Q+ c + Q+ h - Q- d . (4) 

The vertically integrated viscous heating rate is given by 
Qvsc = { T '^)' v y with t as given above. 

As the entropy gradient in A DAFs decrease s outward, 
ADAFs are convectively unstable (Honma 1996). Thus we 
additionally consider a bulk energy transport due to a con- 
vective energy flux P t rb with the associated net heating rate 
given by Q^ lb = — V-ftrb- A natural assumption is to set 
-Ftrb proportional to the local entropy gradient, i.e., 

F trb = - KT STVs, (5) 

with kt an effective di f fusion coefficient. Following llHpnmal 
ll99fJ: IManmoto et alJ l200d: IManmoto fc Katol l2000f) we 

parametrise kt in analogy with the a-viscosity prescrip- 
tion as kt = otTCsH. In general, «t ^ a is expected, 
since convective turbulence contributes to the viscosity, but 
not all sources of viscos ity produce a bulk energy transport 
iRvu fc Goodmanll992l) . The rate at which heat is advected 
with the flow is Q~ dv = ET(v- V)s = V-Ev + P g V-v. The 
radiative cooling rate is Q~ &d = 2cpHKp(aRT 4 — U/2H), 
with up the Planck-mean opacity for free- free absorption. 
Finally, the radiation energy equation is 



d t U + Q dm = Qtc - Qcsc, (6) 

where Qt TC = Q^d- Sinks for the radiative energy density 
originate from radiative losses through the upper and lower 
disk surfa ces QZ, r and from radial d iffusion Q dls ~ V-f^d. 
Following IManmoto fc Katol (12000), we write 

Frad = -t^— VnUen, (7) 

where hr is the Rosseland-mean opacity including electron 
scattering and free- free absorption. The equations are closed 
by an equation of state for the pressure, P = P g + P ra d- We 
apply P g = (7 — 1)E, with 7 the ratio of specific heats, 
together with P ra d = (7/3. 

We discretise the equations on a non-uniform radial 
mesh, and use a time semi- implicit, operator-split scheme to 
advance the solution forward in time. The most restrictive 
constraint on the time-step in explicit schemes is due to ra- 
diative and viscous terms as a result of the associated time- 
scales. We avoid this problem by treating these terms im- 
plicitly. Implicit methods for finite volume discretisation of 
our equations result in a large sparse matrix, where only the 
neighbourhood of the diagonal is populated. The coupling 
between individual equations gives raise to a few off-diagonal 
entries. If this coupling is weak, e.g. if the time increment is 
small, these off-diagonal terms may be neglected at low or- 
der and the resulting matrix can be inverted very efficiently. 
In fact, we first solve each equation independently, then 
update the dynamical variables and finally recompute the 
global defects with the full system. The procedure is now re- 
peated until global con vergence is reached. Th i s sche me was 
originally proposed bv lHuieirat fc Rannacheil il998l) and is 
naturally advantageous for moderately coupled systems. If 
the coupling is too strong, the efficiency in terms of achiev- 
able physical simulation time over computational CPU-time 
is very low, and may even drop below that of an explicit 
scheme. 

In this letter we restrict attention to the model param- 
eters (rhi, a, ot), where rhi is the mass accretion rate of the 
initial model in Eddington units. The black hole mass and 
the adiabatic index are held fixed at M = 10M© and 7 = 1.5 
in all m odels. As initial data we take a self-similar ADAF 
profile of lNaravan fc yJ jl994l) with the advection efficiency 
parameter / a step- function, i.e., with / dropping from unity 
(initial ADAF) to zero (initial SSD) at an arbitrarily cho- 
sen initial transition radius. At the outer boundary, located 
at 2000R g - where R g = GM/c 2 is the gravitational radius 
- we fix the accretion rate to rhi and the specific angular 
momentum to the local Keplerian value. Vanishing gradi- 
ent boundary conditions are applied to the remaining vari- 
ables. At the inner boundary at 2.2 R g , we apply free outflow 
boundary conditions. 



3 A FIDUCIAL SOLUTION 

The transition radius R tT is defined as the location where 
the vertical optical depth is equal to unity. Variables at the 
transition radius are denoted by the subscript 'tr', and the 
subscript 'K' refers to Keplerian motion. Time is measured 
in units of the Keplerian orbital period taken at the initial 
transition radius, P^, unless explicitly noted otherwise. Our 
unit of length is the gravitational radius R g . In Fig.^we plot 
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Figure 1. Radial profile of the contributions to the gas energy 
balance (eq. |1J for our fiducial solution. Note that Qj^ b does act 
as a cooling or heating rate at different radii, i.e. it changes sign 
in the transition region. The dots superimposed on the Q~^ rh curve 
indicate the spatial resolution. 



On top of these rapid oscillations we observe a slow 
outward drift of the transition radius and consequently, the 
SSD shrinks at this very rate while the ADAF grows. For 
the fiducial model introduced above, R tI grows from initially 
55 R g to finally 123 R g at t = 92 Pg- The associated drift 
velocity decreases with increasing R tT or, accordingly, with 
increasing time. A drifting behaviour of this type is per- 
ceived in all our simulations. However, the overall bimodal 
disc structure is never disrupted and the transition region is 
always narrow in the sense defined above. 

At t — 92 Pjj the time-step of our simulation dropped 
drastically to very small values, and the simulation was 
halted. This behaviour was seen in all our simulations 
jGraciall2003) . Interestingly, the values we measure for the 
transition radius at that time is nearly equal to the max- 
imum transition radius predicted for a stationary bimodal 
disc model with the same par ameters and boundary condi- 
tions (Ma nmoto fc Katoll200Gf) . Near this maximum transi- 
tion radius the equations become very stiff, and thus the 
convergence of our numerical method decreases while the 
time-step freezes. The issue of the drift as such is further 
discussed in section l4~2l 



the contributions to the energy balance for a fiducial model 
with parameters (rhi,a, qt) = (0.021,0.4,0.4). The transi- 
tion in the initial conditions was located at R tI = 55R g . 
These parameters were chosen to facilitate comparison with 
other published results. 

We generally identify three domains within the flow ac- 
cording to the magnitude of the different terms in the sta- 
tionary energy equation even though our solutions are 
intrinsically time-dependent and cannot be directly com- 
pared with steady state profiles. With increasing radial dis- 
tance from the central black hole we identify the inner 
ADAF where the optical depth is very low. Next we de- 
fine the transition region as the zone where the rotation is 
super-Keplerian. Super-Keplerian rotati on is present when- 
ever the transition region is narrow (cf. lAbramowicz et ail 
Il99ct) . Between the outer edge of the transition region and 
the outer edge of the disc we observe the SSD where the 
optical depth is very large. The radial profiles are never sta- 
tionary, as discussed below. Nevertheless, the time-averaged 
profiles are rathe r simi lar to the stationary solutions of 
iManmoto fc Katol J2000). 



4 EVOLUTION OF THE TRANSITION 
RADIUS 

The transition radius -Rtr derived from our solutions is a 
function of time. We find that R tT is always located within 
the transition region. After a short initial relaxation phase, 
we generally observe a quasi-periodic time variation of R tI , 
superimposed on a much slower outward directed drift. This 
behaviour is shown in Fig. |5Ja). 

The mass flow rate through the transition region is 
not constant in time, but modulated with the same quasi- 
periodic pattern as R tI . The transition region tends to 
abruptly deposit cold gas into the ADAF. Whenever this 
happens, the transition region grows at the expense of the 
ADAF. 



4.1 The oscillation cycle 

The period of an oscillation cycle At osc increases with time. 
In Fig. Eta) we plot the innermost location and the outer- 
most location which the transition radius takes during the 
complete evolution of our fiducial model. The outer loca- 
tion is treated as the generic transition radius and fitted by 
an analytical function lZt T {t). The innermost location was 
in turn considered as the perturbed transition radius in the 
course of the oscillation cycle and used to derive the oscilla- 
tion period At osc . 

As shown in Fig. IHb), a strong correlation is found 
between At osc and the generic transition radius TZ tI (t). Fur- 
thermore, Atosc is very close to (within a factor of order 
unity) the Keplerian orbital period at the generic transi- 
tion radius, -FkC^tr^)), independent of model parameters 
(rhi, a, err). However, the width of the transition region is a 
function of the adiabatic index 7 (M anmoto fc Katal200fJ) 
and therefore we expect it to modify the proportionality fac- 
tor between Ai osc and Jk(7?.tr(t))- 

A necessary condition for a narrow transition region is 
super-K eplerian rotation within the transition region, i.e. 
I > Ik i Abramowicz et ail 1199 81 . In our calculations the 
specific angular momentum profile decreases outward and 
smoothly matches the Keplerian rotation profile of the SSD. 
However, if I decreases with increasing radius, 8rI < 0, the 
flow is Rayleigh unstable. Consequently, a narrow ADAF- 
SSD transition region is the location of putative instability. 

Outward decreasing specific angula r momen t um p rofiles 
have previously been ob tained by e.g. iHonmal (1996) and 
IManmoto fc Katd feOOOft as seen f rom their figures. The spe- 
cific angular momentum profile of lAbramowicz et alJ lll99St) 
does not decrease outward. They used an adiabatic equa- 
tion of state, which yields broader transition regions, i.e. I 
changes slower near the outer edge of the transition region 
and matches the SSD profile without actually decreasing. 

Local stability analysis of a rotating flow predicts linear 
perturbations to oscillate with frequency uj that is bounded 
by the epicyclic frequency k. Perturbations with uj 2 » 
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Figure 2. (a) temporal evolution of the transition radius Rt T (dotted). The crosses + and stars X indicate the innermost and outermost 
location of the transition radius. The latter is fitted by an analytical function 72-tr(0 (solid), (b) dependence of the oscillation cycle period 
At osc on the transition radius 7^tr for our fiducial model. In both panels time is measured in units of the Keplerian orbital period 
at the initial transition radius located at Rt T = 55 R g . 



will p referably grow or decay (|Kato fc ManmotcfeOOftlKatol 
120011 . for a review). If k 2 is negative the flow is dynamically 
unstable. If k 2 is positive the perturbation wi l l oscil late with 
frequency tv 2 near K^ax- 

iKato fc Manmotol feOOCft further 
argue that the oscillations will be trapped within the transi- 
tion region, since they are reflected at the transition region 
boundaries. 

In our solutions we find an upper limit ^max/^K ~ 2, 
i.e., the oscillation period At OBC /P£ has a lower limit l/\/2 
of order unity, consistent with our numerical results and 
with our proposition that the oscillations have a physical 
nature. On the other hand, we note that K^ax sensitively 
depends on the angular momentum gradient in the transi- 
tion region, which is rather poorly resolved in our current 
models. Therefore, we cannot conclude from our simulations 
that At oac is always bounded by the above limit. 

We mentioned earlier that perturbations generated 
within the transition region are expected to be effectively 
trapped within the same. In reality, this simple picture does 
not survive the confrontation with numerical experiments 
entirely unharmed. Perturbations from within the transition 
region do back-react with the flow outside. Fig.[3]illustrates 
the dynamical evolution of the flow during an oscillation 
cycle at high temporal resolution. 

The oscillation cycle appears to be initiated by a pertur- 
bation as seen, e.g., in the radial momentum flux density, ap- 
proaching the inner boundary of the transition region. The 
perturbation leaks into the ADAF, cold matter with high 
angular momentum detaches from the transition region and 
penetrates the ADAF as seen in Fig.|3] The ejected blob car- 
ries nearly constant angular momentum, such that the re- 
gion of super-Keplerian rotation, i.e., the transition region, 
grows far into the ADAF. When the perturbation reaches 
the sonic radius, it is partly reflected as an outgoing wave 
and sweeps up hot plasma towards the transition. In contin- 



uation the ADAF relaxes slowly back into the unperturbed 
state. This causes high variability not only of the transition 
region, but also of the inner ADAF, where most of the X- 
ray flux is expected to originate. Hence the evolution of our 
models is strongly affected by the existence of oscillations in 
the transition region. 

4.2 Radial drift or the strong ADAF principle 

The outward drift of the transition region which we see in 
our simulations does also require interpretation. It is well 
known that, for a given accretion rate, the radial extent of 
the ADA F (if existent ) is restricted by an upp er limit R T 
(e.g. IChen et alJ Il995t lAbramowicz et al.lfl995h . At every 
radius R > R T , the optical depth exceeds unity and radiative 
cooling is efficient, such that the ADAF solution branch is no 
longer available. For certain R < R T , the ADAF and SSD- 
like solutions are both available in principle. From entropy 
considerations we expect, that nature realizes the ADAF 
whenever possibl e, a proposition refer red to as the strong 
ADAF principle jNaravan fc Yill995rj) . 

Our models do indeed seem to support this hypothe- 
sis. The transition radius 7?.tr(t) does always drift outward, 
thereby increasing the total entropy of the entire solution. 
The drift is stopped once R tI has reached the outermost lo- 
cation, i?trb, that is allowed from stationary analysis. Note 
that J? t rb is smaller than R T because turbulent convective 
energy transport is taken into account in our models. 

For a given set of parameters, the specific entropy of the 
ADAF is, at any radius of interest, larger than the specific 
entropy of the SSD. The particular form applied for turbu- 
lent convective energy transport (namely along the negative 
entropy gradient), implies that the transition region is al- 
ways supplied with energy from the inner ADAF. As long 
as this energy surplus cannot be completely shed by radia- 
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Figure 3. Dynamics of the oscillations in the ADAF-SSD transition region. On each panel radius R = (2R g ■ ■ ■ 2000R g ) increases from 
left to right on a logarithmic scale, while time t = (0 • ■ ■ 0.92 P^ r (62 R g )) runs from bottom to top on a linear scale. The individual panels 
show logarithm of temperature logT (left), and logarithm of vertically integrated radial momentum density log(|St)ij|) (right). T is 
measured in units of (7 — l)c 2 (j,m p /kB and TiVr in units of M Edd/(2vrR g ). The outer black contours indicate the zone of super-Kcplerian 
rotation, i.e. the transition region. The inner one, the transition radius Rt T . 



tive losses, the transition region is heated to nearly the local 
virial temperature and the flow switches to the ADAF mode. 
As a consequence, the SSD is slowly evaporated radially at 
its inner edge. This process does not stop until radiative 
energy shedding in the transition region completely domi- 
nates energy supply by convection - exactly the definition 
of the outermost allowed transition radius Rub- Our models 
do indeed evolve in gene ral within the parameter space of 
iManmoto fc Katd fcOOd. see their Figs. 8, 9 and 10). 



5 CONCLUDING REMARKS 

We have modelled the time evolution of bimodal accretion 
discs where an outer cooling dominated disc matches an in- 
ner ADAF-type flow through a narrow transition region. Our 
models show an intrinsically time-dependent behaviour that 
can be understood from the dynamical instabilities predicted 
for stationary ADAF-SSD transitions. Specifically, oscilla- 
tions with frequency close to Qk are excited in the transi- 
tion region. We associate these oscillations with a mutual 
local instability that originates from the radially decreasing 
specific angular momentum profile, a characteristic feature 
of any narrow-transition bimodal disc. 

We observe that the waves excited in the transition re- 
gion leak into the inner ADAF, but do not essentially prop- 
agate outward into the SSD. These oscillations give rise to 
rapid periodic modulations of the ADAF and the transi- 
tion region. We estimated the total disk luminosity by in- 
tegrating the radiative cooling rate over the disk surface 
and again find a periodic modulation with a variability of a 
few percent. This could possibly explain the QPO s obse rved 
in X-Ray transient systems (e.g. |Nowak Lehdll998L and 
references therein) . iKato fc Manmotol i200d) estimated the 
oscillation frequency in bimodal disc transitions to a small 
fraction ~ 10 -3 of Ok, albeit under the assumption of small 
frequency and using a rather crude estimate for the disc 



structure. Our numerical results show frequencies very close 
to Ok at the transition radius. 

We furthermore observe that the initial transition ra- 
dius does slowly drift to larger radii, i.e., even though a 
stationary transition is allowed for a wide range of radii, the 
disc appears to favour the outermost possible transition ra- 
diu s. This fits well into th e strong ADAF principle hypothe- 
sis ijNaravan fc Yill995bl) . Within our modelling, convective 
energy transport is the key to make this drift possible, i.e., 
to evaporate the SSD radially until the transition reaches 
its outermost allowed location. But even if we neglect con- 
vective energy transfer as modelled in this paper, the large 
entropy gradient across the transition still persists. Signifi- 
cant transfer of heat from the ADAF to the SSD is a nec- 
essary condition for the existence of a stationary transition 
front. Whether turbulence is the relevant physical mecha- 
nism for this heat flux and whether the shape of the transi- 
tion front is actually close to cylindrical cannot be addressed 
by 1-dimensional modeling. The resolution of this issue does 
require at least 2-dimensional simulations in the parame- 
ter range that allows for a transition in the 1-dimensional 
stationary case. This region of the parameter space, how- 
ever, has not been addressed in a 2-dimensional simulation 
by Hujeirat & Camenzind (20 00). The gap between ver- 
tical disc evaporation m odels iMever fc Mever-Hofmeisteil 
11994 iMever et al]l200Cft and the present 1-dimensional ra- 
dial transition physics therefore remains to be closed in a 
future work. 
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